################################################################################
#################   MAIN ARTICLE: MODELS, TABLES AND FIGURES  ##################
################################################################################

cat("\014") 
library(car)
library(dplyr)
library(lme4)
library(texreg)
library(ggplot2)
library(scales)
library(sjPlot)
library(lattice)
library(forcats)
library(gridExtra)
library(ggpubr)
library(effects)
library(ggeffects)
library(gridExtra)
library(labelled)
library(ggplotify)
options(scipen=999)

rm(list=ls())
load("lapop_pela_select_cent_wVH.RData"); baseF=base2; rm(base2)


###########--------------------- MAIN MODELS  ----------------------############

#Centering strategy based on: Ahnalee M. Brincks, Craig K. Enders, Maria M. Llabre, 
#Rebecca J. BulotskyShearer, Guillermo Prado & Daniel J. Feaster (2017) Centering 
#Predictor Variables in Three-Level Contextual Models, Multivariate Behavioral 
#Research, 52:2, 149-163, DOI: 10.1080/00273171.2016.1256753


###### Within survey and contextual country models (CGM1/CGM2) (Table 2) #######
cgm0=lmer(conf_parlamnt ~ 
             eval_econpais_Dummy_gm+ izq_gm+ voto_presid_gm+ conf_interpers_gm+
             eduniv_gm+ urbano_gm+ mujer_gm+ edad_gm+ casado_gm + 
             wave+ (1|pais.ano)+(1|paisF), data=baseF)

cgm1=update(cgm0, . ~ . +   fampol_gm+ fampol_prom)
cgm2=update(cgm0, . ~ . +   reeleccion_gm+ reeleccion_prom)
cgm3=update(cgm0, . ~ . +   fampol_gm+ fampol_prom+ reeleccion_gm+ reeleccion_prom)
screenreg(list(cgm1,cgm2,cgm3), stars=c(0.1,0.05,0.01))

#Controlling by contextual variables (V-DEM and WB)  
cgm1_a = update(cgm1, . ~ . + e_gdppc_gm + v2x_corr_gm + desempleo_gm)
cgm2_a = update(cgm2, . ~ . + e_gdppc_gm + v2x_corr_gm + desempleo_gm)
cgm3_a = update(cgm3, . ~ . + e_gdppc_gm + v2x_corr_gm + desempleo_gm)
screenreg(list(cgm1_a,cgm2_a,cgm3_a), stars=c(0.1,0.05,0.01))





############  Split-sample models (conditional effects)  #############

############ Interaction models (conditional effects) (Table 3)  ###############

quantile(baseF$fampol_prom, na.rm=T)
baseF$fampol_promF <- cut(baseF$fampol_prom, include.lowest=TRUE, 
                          breaks=quantile(baseF$fampol_prom, c(0,.5,1), na.rm=T))
levels(baseF$fampol_promF) <-c("Low", "High"); table(baseF$fampol_promF)

### Reelection (Low and High)
quantile(baseF$reeleccion_prom, na.rm=T)
baseF$reeleccion_promF<-cut(baseF$reeleccion_prom, include.lowest=TRUE,
                            breaks=quantile(baseF$reeleccion_prom, c(0,.5,1), 
                                            na.rm=T))
levels(baseF$reeleccion_promF) <-c("Low", "High")
table(baseF$reeleccion_promF)

#With interactions 
cgm1_ai=lmer(conf_parlamnt ~ 
               eval_econpais_Dummy_gm+ izq_gm+ voto_presid_gm+ conf_interpers_gm+
               eduniv_gm+ urbano_gm+ mujer_gm+ edad_gm+ casado_gm + 
               wave+ e_gdppc_gm + v2x_corr_gm + desempleo_gm + 
               fampol_gc*fampol_promF +  
               (1|pais.ano) + (fampol_gc|paisF), 
             data=baseF)

cgm2_ai=lmer(conf_parlamnt ~ 
               eval_econpais_Dummy_gm+ izq_gm+ voto_presid_gm+ conf_interpers_gm+
               eduniv_gm+ urbano_gm+ mujer_gm+ edad_gm+ casado_gm + 
               wave+ e_gdppc_gm + v2x_corr_gm + desempleo_gm +
               reeleccion_gc*reeleccion_promF +  
               (1|pais.ano) + (reeleccion_gc|paisF), 
             data=baseF)


#############################      TABLES      #################################
library(vtable)
library(dplyr)
library(kableExtra)
library(fBasics)


################------------ Table 1: Descriptives -------------################

tab1=dplyr::select(baseF, conf_parlamnt,#dependent variable
                   familiar_pol_media_w, reeleccion,#independent variables level 2
                   fampol_prom, reeleccion_prom,#independent variables level 3
                   e_gdppc, v2x_corr, desempleo,#controls
                   eval_econpais_Dummy, conf_interpers, voto_presid, izq,#controls
                   eduniv, urbano, mujer, edad, casado)#controls 

###
tab1=as.data.frame(t(basicStats(tab1)))
tab1=tab1[,-c(5,6,9:13,15:16)]
tab1=tab1[,c(3:7,1,2)]
tab1[,1:5] <- round(tab1[,1:5], digits = 2) 
rownames(tab1)<-c("Trust in the Parlament", 
                  "Relatives in Politics", "Reelection Rate",
                  "Relatives in Politics ", "Reelection Rate ",
                  "GDP Percapita ", "Corruption Index", "Unemployment Rate",
                  "Perception of National Economy", "Social Trust", 
                  "Turnout on Presidential Election", "Left Ideology",
                  "Tertiary Education", "Urban", "Female", "Age", "Married")

kbl(tab1, booktabs = T,
    col.names = c("Min", "Max","Mean","Median","SD", "Obs", "NA"),
    caption = "EstadC-sticos descriptivos de variables") %>% 
  kable_styling(latex_options = c("center", "basic"), font_size = 12)%>%
  pack_rows("Dependent variable", 1,1) %>%
  pack_rows("Independent variables (survey)", 2, 3) %>%
  pack_rows("Independent variables (country)", 4, 5) %>%
  pack_rows("Control variables (survey)", 6, 8) %>%
  pack_rows("Control variables (individual)", 9, 17)



#################------------ Table 2: Main models -------------################

htmlreg(list(cgm1_a,cgm2_a,cgm3_a), stars=c(0.1,0.05,0.01),
        single.row = F, caption.above = T, digits= 2,
        file = "tabla2.doc",
        caption = "Elite Social Closure Effects on Trust in the Parliament at the survey and country levels.",
        custom.model.names = c("1","2","3"),
        custom.coef.map= list("(Intercept)"= "Intercept", 
                              "eval_econpais_Dummy_gm"="Perception of the economy(cgm)", 
                              "izq_gm"="Left ideology (cgm)",
                              "voto_presid_gm"="Turnout on Presidential Election (cgm)", 
                              "conf_interpers_gm"="Social Trust (cgm)", 
                              "eduniv_gm" = "Terciary education (cgm)",
                              "urbano_gm" = "Urban (cgm)",
                              "mujer_gm" = "Female (cgm)",
                              "edad_gm" = "Age (cgm)",
                              "casado_gm" = "Married (cgm)",
                              "e_gdppc_gm"="GDP Percapita (cgm)", 
                              "v2x_corr_gm"="Political Corruption (cgm)",
                              "desempleo_gm"="Unemployment Rate (cgm)",
                              "fampol_gm"="Proportion of Relatives in Politics(cgm)", 
                              "reeleccion_gm"="Reelection Rate (cgm)",
                              "fampol_prom"="Proportion of Relatives in Politics",
                              "reeleccion_prom"="Reelection Rate"),
        groups = list("Control variables (individual-level)"=2:10, 
                      "Control variables (survey-level)"= 11:13, 
                      "Social closure (survey-level)"=14:15,
                      "Social closure (country-level)"=16:17),
        include.variance = FALSE, 
        custom.note = "%stars. cgm = centered on the grand mean. Models include a linear parameter to capture the time trend. Data from LAPOP, PELA-USAL V-DEM and ILO.")


###########------------ Table 3: Models with interactions ------------###########
screenreg(list(cgm1_ai,cgm2_ai), stars=c(0.1,0.05,0.01))
htmlreg(list(cgm1_ai,cgm2_ai), stars=c(0.1,0.05,0.01),
        single.row = F, caption.above = T, digits= 2,
        file = "tabla3.doc",
        caption = "Elite Social Closure Effects on Trust in the Parliament at the Survey Level Moderated for the Country Level of Closure.",
        custom.model.names = c("4","5"),
        custom.coef.map= list("(Intercept)"= "Intercept", 
                              "eval_econpais_Dummy_gm"="Perception of the economy(cgm)", 
                              "izq_gm"="Left ideology (cgm)",
                              "voto_presid_gm"="Turnout on Presidential Election (cgm)", 
                              "conf_interpers_gm"="Social Trust (cgm)", 
                              "eduniv_gm" = "Terciary education (cgm)",
                              "urbano_gm" = "Urban (cgm)",
                              "mujer_gm" = "Female (cgm)",
                              "edad_gm" = "Age (cgm)",
                              "casado_gm" = "Married (cgm)",
                              "e_gdppc_gm"="GDP Percapita (cgm)", 
                              "v2x_corr_gm"="Political Corruption (cgm)",
                              "desempleo_gm"="Unemployment Rate (cgm)",
                              "fampol_gc"="Proportion of Relatives in Politics (cwc)", 
                              "fampol_promFHigh"="High Relatives in Politics",
                              "fampol_gc:fampol_promFHigh" = "Interaction",
                              "reeleccion_gc"="Re-election Rate (cwc)",
                              "reeleccion_promFHigh" = "High Re-election Rate",
                              "reeleccion_gc:reeleccion_promFHigh"= "Interaction"),
        groups = list("Control variables (individual-level)"=2:10, 
                      "Control variables (survey-level)"= 11:13, 
                      "Elite Social closure"=14:19),
        include.variance = FALSE, 
        custom.note = "%stars. cgm = centered on the grand mean, cwc= centered within context, Models also include country fixed-effects and a linear parameter to control the time trend. These estimates were omitted from the table. Data from LAPOP, PELA-USAL V-DEM and ILO.")


#############################      FIGURES      ################################


############------- Figure 1: Political Trust in Latin America ------###########

conf.pais.ano <- baseF %>% group_by(paisR, wave) %>% 
  summarise(confpol = mean(conf_parlamnt, na.rm=T))
conf.pais.ano <- conf.pais.ano %>% group_by(paisR) %>% 
  mutate(mean.pais = mean(confpol, na.rm = T))
conf.pais.ano <- conf.pais.ano %>% group_by(wave) %>% 
  mutate(mean.wave = mean(confpol, na.rm = T))
conf.pais.ano$mean.general <- mean(conf.pais.ano$confpol, na.rm=T)
conf.pais.ano$paisR <-car::recode(conf.pais.ano$paisR, "'Brasil'='Brazil';
                                  'M??xico'='Mexico'; 'Panam??'='Panama';
                                  'Per??'='Peru'; 'Panam??'='Panama'; 
                                  'Rep. Dominicana'='Dominican Rep.'")
set_theme(theme_bw())
ggplot(data=conf.pais.ano, aes(x=wave, y=confpol))+
  geom_line()+geom_point(size=1) +
  geom_hline(aes(yintercept = mean.pais, col='red'), linetype='dashed') +
  facet_wrap(~paisR, ncol=3)+
  scale_y_continuous(name="", 
                     limits=c(1.9,5), breaks=c(2,3,4,5))+
  scale_x_continuous(name="", 
                     breaks=c(2004,2006,2008,2010,2012,2014,2016,2018))+
  theme(axis.title.x = element_blank(),
        axis.text.y = element_text(size=7.5),
        axis.text.x = element_text(size=7),
        legend.position = 'none')



########---- Figure 2 and 3: Elite Social Closure in Latin America -----########

horiz <- baseF %>% group_by(paisR, wave) %>% 
  summarise(fampol = mean(familiar_pol_media_w, na.rm=T),
            reeleccion=mean(reeleccion, na.rm=T))

means <-baseF %>% group_by(paisR) %>% 
  summarise(fampol_gm = mean(familiar_pol_media_w, na.rm=T),
            reeleccion_gm=mean(reeleccion, na.rm=T))

horiz <-merge(horiz, means, by="paisR")
horiz$paisR <-car::recode(horiz$paisR, "'Brasil'='Brazil';
                                  'M??xico'='Mexico'; 'Panam??'='Panama';
                                  'Per??'='Peru'; 'Panam??'='Panama'; 
                                  'Rep. Dominicana'='Dominican Rep.'")

#Figure 2: Relatives in politics
set_theme(theme_bw())
ggplot(data=horiz, aes(x=wave, y=fampol))+ geom_line()+geom_point(size=1)+
  geom_hline(data =horiz, aes(yintercept=fampol_gm, color= "red"), 
             linetype="dashed")+
  facet_wrap(~paisR, ncol=3)+
  scale_y_continuous(name="", limits=c(0.2,0.8), 
                     breaks=c(0.3,0.5,0.7))+
  scale_x_continuous(name="", 
                     breaks=c(2004,2006,2008,2010,2012,2014,2016,2018))+
  theme(axis.title.x = element_blank(),
        axis.text.y = element_text(size=7.5),
        axis.text.x = element_text(size=7),
        legend.position = 'none')


#Figure 3: Reelection Rate
set_theme(theme_bw())
ggplot(data=horiz, aes(x=wave, y=reeleccion))+ geom_line()+geom_point(size=1)+
  geom_hline(data = horiz, aes(yintercept=reeleccion_gm, color = "red"), 
             linetype="dashed")+
  facet_wrap(~paisR, ncol=3)+
  scale_y_continuous(name="", limits=c(0,0.9), 
                     breaks=c(0.2,0.5,0.8))+
  scale_x_continuous(name="", 
                     breaks=c(2004,2006,2008,2010,2012,2014,2016,2018))+
  theme(axis.title.x = element_blank(),
        axis.text.y = element_text(size=7.5),
        axis.text.x = element_text(size=7),
        legend.position = 'none')



######--- Figure 4: Latin American countries by elites' social closure ---######

baseF$paisF <-car::recode(baseF$paisF, "'Brasil'='Brazil'; 'M??xico'='Mexico'; 
                          'Panam??'='Panama'; 'Per??'='Peru'; 'Panam??'='Panama'; 
                          'Rep. Dominicana'='Dominican Rep.'")

rp <- baseF %>% group_by(paisF) %>% summarise(fampol=mean(fampol_prom))
rp$fampolF <-cut(rp$fampol, include.lowest=TRUE,
                 breaks=quantile(rp$fampol, c(0,.5,1), na.rm=T))
levels(rp$fampolF) <-c("Low", "High")
table(rp$fampolF)

rr <- baseF %>% group_by(paisF) %>% summarise(reeleccion=mean(reeleccion_prom))
rr$reeleccionF <-cut(rr$reeleccion, include.lowest=TRUE,
                     breaks=quantile(rr$reeleccion, c(0,.5,1), na.rm=T))
levels(rr$reeleccionF) <-c("Low", "High")
table(rr$reeleccionF)

#Box plot
set_theme(theme_bw())
rpF <-ggplot(rp, aes(x = "", y = fampol)) + 
  scale_colour_manual("fampol", values=c("grey41", "black"))+
  labs(y = "Relatives in Politics", x = "") + labs(linetype = "")+
  geom_jitter(aes(color=fampolF), position = position_jitter(seed = 1),size=1.1)+
  ggrepel::geom_text_repel(aes(label=paisF, color=fampolF),
                           position = position_jitter(seed=1), size=3.2)+
  geom_hline(yintercept = median(rp$fampol), linetype="dashed", 
             color="red", size=0.7)+
  theme(legend.position = "none",
        legend.title = element_blank(),
        legend.text = element_text(size=9),
        axis.text.x = element_text(size=8),
        axis.title.x = element_text(size=10),
        axis.ticks.y = element_blank())+
  coord_flip()
rpF

set_theme(theme_bw())
rrF <-ggplot(rr, aes(x = "", y = reeleccion)) + 
  scale_colour_manual("reeleccion", values=c("grey41", "black"))+
  labs(y = "Reelection Rate", x = "") + labs(linetype = "")+
  geom_jitter(aes(color=reeleccionF), position = position_jitter(seed =1),size=1.1)+
  ggrepel::geom_text_repel(aes(label=paisF, color=reeleccionF), 
                           position = position_jitter(seed = 1),size=3.2)+
  geom_hline(yintercept = median(rr$reeleccion), linetype="dashed", 
             color="red", size=0.7)+
  theme(legend.position = "none",
        legend.title = element_blank(),
        legend.text = element_text(size=9),
        axis.text.x = element_text(size=8),
        axis.title.x = element_text(size=10),
        axis.ticks.y = element_blank())+
  coord_flip()
rrF

ggarrange(rpF,rrF, ncol =1, nrow =2 ,common.legend =TRUE, legend="bottom")



###########------------------ Figure 5: Interactions ----------------###########

#coefficients
fampol_l <- fixef(cgm1_ai)[15]
fampol_h <- fixef(cgm1_ai)[15] + fixef(cgm1_ai)[17]
reelec_l <- fixef(cgm2_ai)[15]
reelec_h <- fixef(cgm2_ai)[15] + fixef(cgm2_ai)[17]
coefs<-as.numeric(c(reelec_l,reelec_h,fampol_l,fampol_h))

#standard errors
sefampol_l<-sqrt(vcov(cgm1_ai)[15, 15])
sefampol_h<-sqrt(vcov(cgm1_ai)[15, 15] + vcov(cgm1_ai)[17, 17]+ 2*vcov(cgm1_ai)[15, 17])
sereelec_l<-sqrt(vcov(cgm2_ai)[15, 15])
sereelec_h<-sqrt(vcov(cgm2_ai)[15, 15] + vcov(cgm2_ai)[17, 17]+ 2*vcov(cgm2_ai)[15, 17])
ses<-c(sefampol_l, sefampol_h, sereelec_l, sereelec_h)

#CI 95% (df= 97)
li95_fampoll<-fampol_l - (1.98*sefampol_l)
ls95_fampoll<-fampol_l + (1.98*sefampol_l)
li95_fampolh<-fampol_h - (1.98*sefampol_h)
ls95_fampolh<-fampol_h + (1.98*sefampol_h)

#CI 90% (df= 97)
li90_fampoll<-fampol_l - (1.67*sefampol_l)
ls90_fampoll<-fampol_l + (1.67*sefampol_l)
li90_fampolh<-fampol_h - (1.67*sefampol_h)
ls90_fampolh<-fampol_h + (1.67*sefampol_h)

#CI 95% (df= 97)
li95_reelecl<-reelec_l - (1.98*sereelec_l)
ls95_reelecl<-reelec_l + (1.98*sereelec_l)
li95_reelech<-reelec_h - (1.98*sereelec_h)
ls95_reelech<-reelec_h + (1.98*sereelec_h)

#CI 90% (df= 97)
li90_reelecl<-reelec_l - (1.67*sereelec_l)
ls90_reelecl<-reelec_l + (1.67*sereelec_l)
li90_reelech<-reelec_h - (1.67*sereelec_h)
ls90_reelech<-reelec_h + (1.67*sereelec_h)

var<-c("Re-election rate", "Re-election rate","Relatives in politics", 
       "Relatives in politics")
group<-c("Low","High","Low","High")
liminf90<-c(li90_reelecl, li90_reelech, li90_fampoll, li90_fampolh)
liminf95<-c(li95_reelecl, li95_reelech, li95_fampoll, li95_fampolh)
limsup90<-c(ls90_reelecl, ls90_reelech, ls90_fampoll, ls90_fampolh)
limsup95<-c(ls95_reelecl, ls95_reelech, ls95_fampoll, ls95_fampolh)
df4<-as.data.frame(cbind(var,group,coefs,liminf90,liminf95,
                         limsup90,limsup95))
df4$coefs<-as.numeric(df4$coefs)
df4$liminf90<-as.numeric(df4$liminf90);df4$liminf95<-as.numeric(df4$liminf95)
df4$limsup90<-as.numeric(df4$limsup90);df4$limsup95<-as.numeric(df4$limsup95)

#Figure 5
theme_set(theme_bw())
ggplot(df4, aes(coefs, group)) + geom_point(size=2) +
  geom_linerange(aes(xmin = liminf90, xmax = limsup90), linewidth = 1) +
  geom_linerange(aes(xmin = liminf95, xmax = limsup95), linewidth = 0.5)+
  geom_vline(xintercept = 0, linewidth = 0.5, colour = "red", linetype = "dashed")+
  facet_wrap(~var, nc=1)+
  labs(x="Point Estimate with 90% and 95% Confidence Intervals", y="")+
  theme(axis.title.x = element_text(size=9),
        strip.text.x = element_text(size = 11))